Author

Emi Carlevaro

Published

January 19, 2026

INIT

Code
#required packages to make a package
pacman::p_load(devtools, roxygen2, testthat)
#dependecies
#pacman::p_load(tidyverse, data.table, DT, plotly, gridExtra, latex2exp, here, logr, future.apply, furrr)
#pacman::p_load(dplyr,parallel,sandwich,here,future,logr,profvis,lobstr,parallel)
#https://ourcodingclub.github.io/tutorials/writing-r-package/
#create_package
#https://combine-australia.github.io/r-pkg-dev/functions.html
#load_all()
#run these two lines and then library(GridMaker) in normal rstudio to install package
# load_all loads a package. It roughly simulates what happens when a package is installed and loaded with library().
devtools::load_all()

devtools::install()
Code
library(vars)
Loading required package: MASS
Loading required package: strucchange
Loading required package: zoo

Attaching package: 'zoo'
The following objects are masked from 'package:base':

    as.Date, as.Date.numeric
Loading required package: sandwich
Loading required package: urca
Loading required package: lmtest
Code
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ stringr::boundary() masks strucchange::boundary()
✖ dplyr::filter()     masks stats::filter()
✖ dplyr::lag()        masks stats::lag()
✖ dplyr::select()     masks MASS::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code
library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:MASS':

    select

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout
Code
library(iStabi)
Warning: replacing previous import 'purrr::transpose' by
'data.table::transpose' when loading 'iStabi'
Warning: replacing previous import 'data.table::first' by 'dplyr::first' when
loading 'iStabi'
Warning: replacing previous import 'data.table::last' by 'dplyr::last' when
loading 'iStabi'
Warning: replacing previous import 'data.table::between' by 'dplyr::between'
when loading 'iStabi'
Code
sam <- iStabi::data_sim_biVAR

THE MODEL

True parameters are

\[ \begin{align} y_{1,t} &= \textcolor{blue}{\beta} \,\, y_{2,t} + 0.5 y_{1, t-1} + 0.2 y_{2, t-1} + w_{1, t} \\ y_{2,t} &= \textcolor{blue}{\alpha} \,\, y_{1,t} + 0.1 y_{1, t-1} + 0.4 y_{2, t-1} + w_{2, t} \end{align} \]

\[ \begin{bmatrix}1 & -0.15 \\ 0.75 & 1\end{bmatrix} y_t = \begin{bmatrix}0.5 & 0.2 \\ 0.1 & 0.4\end{bmatrix} y_{t-1} + w_t \]

with \(\beta = 0.15\) and \(\alpha = -0.75\).

ESTIMATION VAR

Code
varModel <- VAR(cbind('i' = sam$y1, 's' = sam$y2), p = 1, type = "const")
Res = residuals(varModel)
Code
# Use for computing the moment function
Omega <- cbind('ii' = Res[, 'i']*Res[, 'i'], 'is' = Res[, 'i']*Res[, 's'],
               'ss' = Res[, 's']*Res[, 's'])

Y <- cbind( 'omega_ii' = Omega[, 'ii'], 'omega_ss' = Omega[,'ss'], 'omega_is' = Omega[,'is'] )

# Compute the moment function at each t. 
# INPUT: a named vector with values for the structural parameters
# OUTPUT: a TxG matrix where each row corresponds to the value of the moment function on that obseration (date)
# an each column the value fot ehfat moment equation
GET_MOM_VALUES <- function(t0, ...) {
  list2env(as.list(t0),envir = environment())
  
  #Y <- cbind( 'omega_ii' = Res[,'i']^2, 'omega_ss' = Res[,'s']^2, 'omega_is' = Res[,'i'] * Res[,'s'] )
  
  #Y %*% c(-theta0['alpha'],
  #        -theta0['beta'],
  #        1 + theta0['alpha'] * theta0['beta'])
  
  #fT 
  as.matrix(-alpha*Omega[, 'ii'] - beta*Omega[, 'ss'] + (1+alpha*beta) * (Omega[, 'is']))
  # fT - lambda
  
  #matrix(fT - mean(fT),
  #       ncol = 1) # (T-varLeng) X 1
  
} 
Code
GET_B <- function(theta0) {
  c(-theta0['alpha'],
    -theta0['beta'],
    1 + theta0['alpha'] * theta0['beta'])
}

BUILD SETS

Critical values

Critical values for \(qLL-\tilde{S}\) for 1 moment condition from table VII Suplementary material, p3 Magnusson Mavroeidis (2014).

THe S-test follows a Chi-Square distribution with 1 degrees of freedom (equatl to the number of moment conditions).

Grid has 1600 rows.

The number of fixed identified parameters is 1.

The number of estimated parameters under the null is 0.

Code
NSIPARAMS = 0 # not needed in the new version of the package
PCONFS = c(0.85, 0.90)
get_CVs <- function() {
  
  CVs <- list('qLLStab'=NA, 'S'=NA, 'genS_qLL'=NA)

  CVs$qLLStab <- filter(iStabi::data_CVs$qLLStab, DG_F == KZ & P_CONF %in% PCONFS)
  CVs$genS_qLL <- filter(iStabi::data_CVs$genS_qLL, 
                         MOMENT_CONDITIONS == KZ & STRONGLY_IDENTIFIED_PARAMETERS == NSIPARAMS & P_CONF %in% PCONFS)
  CVs$S <- tibble('P_CONF' = PCONFS,
                 'CUTOFF_TOP' = map_dbl(PCONFS, ~qchisq(.x, KZMKX)))
  
  CVs
}
#CV$S <- tibble('SIG' = c(0.1, 0.05, 0.01, 0),
#               'VALUE' = c(0, qchisq(.9, KZMKX), 
#                           qchisq(.95, KZMKX), 
#                           qchisq(.99, KZMKX)))

CVs <- get_CVs()

Critical values are:

Code
CVs$S
# A tibble: 2 × 2
  P_CONF CUTOFF_TOP
   <dbl>      <dbl>
1   0.85       2.07
2   0.9        2.71
Code
CVs$qLLStab
# A tibble: 2 × 3
   DG_F P_CONF CUTOFF_TOP
  <dbl>  <dbl>      <dbl>
1     1   0.85       6.46
2     1   0.9        7.17
Code
CVs$genS_qLL
# A tibble: 2 × 4
  MOMENT_CONDITIONS STRONGLY_IDENTIFIED_PARAMETERS P_CONF CUTOFF_TOP
              <dbl>                          <dbl>  <dbl>      <dbl>
1                 1                              0   0.85       7.74
2                 1                              0   0.9        8.59

Data for plotting

Code
build_set <- function(setName) {
  # setName = 'qLLStab'
  thisCV <- CVs[[setName]]
  statMax = filter(thisCV, P_CONF == PCONFS[2])$CUTOFF_TOP
  thisPconf = thisCV$P_CONF
  thisCV = thisCV$CUTOFF_TOP
  # 546 rows
  
  theseCols <- c(PARAMS_CFG$name, setName)
  aSet = Grid[get(setName) <= statMax, ..theseCols][, 
                                                    sigAt := thisPconf[ 1+findInterval(get(setName), thisCV,left.open=TRUE)]]
  aSet
  
}
Code
Sets <- list('qLLStab' = build_set('qLLStab'),
             'S' = build_set('S'),
             'genS_qLL' = build_set('genS_qLL'))

PLOTTING

Grid for plotting

Force the starting values to be a multiple of the step size. This avoids problems and is clear what the true parameter space is.

Plots

True parameter vector $ =(= -0.75, = 0.15)$